Practice Dec 12, 2016
Problem 1
Remember the tomato data set generated by Pepe; we first looked at this when we were working on Chapter 9. 35 accessions for seven species were grown in sun and shade.
Assess whether there is evidence for total length (“totleng”) response to shade and whether this response varies by species. Consider whether including accession (“acs”) or shelf (“shelf”) using adaptive priors improves the model fit.
get the data and format
tomato <- read.csv("Assignment_Chapter_09/TomatoR2CSHL.csv")
names(tomato)## [1] "shelf" "flat" "col" "row" "acs" "trt"
## [7] "days" "date" "hyp" "int1" "int2" "int3"
## [13] "int4" "intleng" "totleng" "petleng" "leafleng" "leafwid"
## [19] "leafnum" "ndvi" "lat" "lon" "alt" "species"
## [25] "who"tomato$species_id <- coerce_index(tomato$species)
tomato$shade <- ifelse(tomato$trt=="L",1,0)
tomato$shelf_id <- coerce_index(tomato$shelf)
tomato$species_id <- coerce_index(tomato$species)
tomato$acs_id <- coerce_index(tomato$acs)
tomato.small <- tomato[,c("shelf_id","acs_id","shade","totleng","species_id")]
summary(tomato.small)## shelf_id acs_id shade totleng
## Min. :1.000 Min. : 1.0 Min. :0.0000 Min. : 13.59
## 1st Qu.:2.000 1st Qu.: 9.0 1st Qu.:0.0000 1st Qu.: 39.25
## Median :3.000 Median :18.0 Median :1.0000 Median : 50.98
## Mean :3.512 Mean :18.1 Mean :0.5089 Mean : 53.70
## 3rd Qu.:5.000 3rd Qu.:27.0 3rd Qu.:1.0000 3rd Qu.: 64.76
## Max. :6.000 Max. :36.0 Max. :1.0000 Max. :129.43
## species_id
## Min. :1.000
## 1st Qu.:2.000
## Median :3.000
## Mean :2.927
## 3rd Qu.:4.000
## Max. :5.000shade effect
mean(tomato.small$totleng[tomato.small$shade==0])## [1] 42.0596m.shade <- map2stan(
alist(
totleng ~ dnorm(mu,sigma),
mu <- alpha +
beta_shade*shade,
alpha ~ dnorm(42,10),
beta_shade ~ dnorm(0,20),
sigma ~ dcauchy(0,1)
),
data = tomato.small,
chains=4
)##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 0.270754 seconds (Warm-up)
## 0.206632 seconds (Sampling)
## 0.477386 seconds (Total)
##
##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 2).
##
## Chain 2, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 0.31679 seconds (Warm-up)
## 0.218748 seconds (Sampling)
## 0.535538 seconds (Total)
##
##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 3).
##
## Chain 3, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 0.305597 seconds (Warm-up)
## 0.217553 seconds (Sampling)
## 0.52315 seconds (Total)
##
##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 4).
##
## Chain 4, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 0.30571 seconds (Warm-up)
## 0.206341 seconds (Sampling)
## 0.512051 seconds (Total)## The following numerical problems occured the indicated number of times after warmup on chain 4## count
## Exception thrown at line 19: normal_log: Scale parameter is 0, but must be > 0! 2## When a numerical problem occurs, the Hamiltonian proposal gets rejected.## See http://mc-stan.org/misc/warnings.html#exception-hamiltonian-proposal-rejected## If the number in the 'count' column is small, do not ask about this message on stan-users.##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 3e-06 seconds (Warm-up)
## 0.000116 seconds (Sampling)
## 0.000119 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]plot(m.shade)
pairs(m.shade)
precis(m.shade)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha 42.09 0.68 40.97 43.15 1745 1
## beta_shade 22.79 0.95 21.30 24.27 1703 1
## sigma 16.03 0.34 15.47 16.54 2867 1shade and species
m.shade.species <- map2stan(
alist(
totleng ~ dnorm(mu,sigma),
mu <- alpha +
beta_sp[species_id] +
beta_shade*shade,
alpha ~ dnorm(42,10),
beta_sp[species_id] ~ dnorm(0,10),
beta_shade ~ dnorm(0,20),
sigma ~ dcauchy(0,1)
),
data = tomato.small,
chains=4
)##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 2.30954 seconds (Warm-up)
## 1.41873 seconds (Sampling)
## 3.72827 seconds (Total)
##
##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 2).
##
## Chain 2, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 2.07273 seconds (Warm-up)
## 1.5323 seconds (Sampling)
## 3.60503 seconds (Total)
##
##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 3).
##
## Chain 3, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 1.97243 seconds (Warm-up)
## 1.20629 seconds (Sampling)
## 3.17872 seconds (Total)
##
##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 4).
##
## Chain 4, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4, Iteration: 2000 / 2000 [100%] (Sampling)
## Elapsed Time: 1.87549 seconds (Warm-up)
## 1.32018 seconds (Sampling)
## 3.19567 seconds (Total)
##
##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 3e-06 seconds (Warm-up)
## 0.00016 seconds (Sampling)
## 0.000163 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]plot(m.shade.species)
pairs(m.shade.species)
precis(m.shade.species,depth = 2)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha 41.11 4.13 34.84 48.10 579 1.01
## beta_sp[1] -1.24 4.22 -7.69 5.84 592 1.01
## beta_sp[2] 0.06 4.22 -6.47 7.00 595 1.01
## beta_sp[3] 2.39 4.22 -4.69 8.80 606 1.01
## beta_sp[4] -9.44 4.28 -15.97 -2.61 609 1.01
## beta_sp[5] 8.35 4.20 2.12 15.59 592 1.01
## beta_shade 23.05 0.94 21.57 24.54 1752 1.00
## sigma 15.16 0.32 14.66 15.68 1953 1.00compare(m.shade,m.shade.species) ## WAIC pWAIC dWAIC weight SE dSE
## m.shade.species 8348.2 6.9 0.0 1 53.08 NA
## m.shade 8458.3 3.1 110.1 0 52.60 21.98levels(tomato$species)## [1] "S. chilense" "S. chmielewskii" "S. habrochaites" "S. pennellii"
## [5] "S. peruvianum"Good evidence for species being important
consider species X shade interaction
m.shade.species.int <- map2stan(
alist(
totleng ~ dnorm(mu,sigma),
mu <- alpha +
beta_sp[species_id] +
beta_shade*shade +
beta_sp_shade[species_id]*shade,
alpha ~ dnorm(53,20),
beta_sp[species_id] ~ dnorm(0,10),
beta_shade ~ dnorm(0,20),
beta_sp_shade[species_id] ~ dnorm(0,10),
sigma ~ dcauchy(0,1)
),
data = tomato.small,
chains=4,
cores=2
)##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 4e-06 seconds (Warm-up)
## 0.000239 seconds (Sampling)
## 0.000243 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]plot(m.shade.species.int)
par(mfrow=c(1,1))
precis(m.shade.species.int,depth=2)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha 41.62 4.23 35.31 48.76 952 1
## beta_sp[1] -4.31 4.37 -11.24 2.52 1013 1
## beta_sp[2] 2.62 4.41 -4.63 9.33 990 1
## beta_sp[3] 3.10 4.35 -3.81 9.92 1003 1
## beta_sp[4] -11.82 4.43 -19.30 -5.14 1041 1
## beta_sp[5] 7.05 4.37 0.12 13.89 983 1
## beta_shade 22.56 4.36 15.16 29.17 904 1
## beta_sp_shade[1] 5.51 4.67 -1.52 13.31 975 1
## beta_sp_shade[2] -5.41 4.61 -12.61 2.20 974 1
## beta_sp_shade[3] -2.04 4.61 -8.73 5.95 1033 1
## beta_sp_shade[4] 4.01 4.71 -3.48 11.28 1075 1
## beta_sp_shade[5] 2.06 4.63 -5.68 9.10 1032 1
## sigma 15.03 0.33 14.48 15.56 2693 1levels(tomato$species)## [1] "S. chilense" "S. chmielewskii" "S. habrochaites" "S. pennellii"
## [5] "S. peruvianum"compare(m.shade,m.shade.species,m.shade.species.int)## WAIC pWAIC dWAIC weight SE dSE
## m.shade.species.int 8335.2 11.0 0 1 53.38 NA
## m.shade.species 8348.2 6.9 13 0 53.08 8.02
## m.shade 8458.3 3.1 123 0 52.60 22.15plot(coeftab(m.shade,m.shade.species,m.shade.species.int))
Reasonable evidence for the interaction, although the 95% posterior intervals overlap among all of the species_shade interactions
Add adaptive prior for accessions
m.shade.species.int.acs <- map2stan(
alist(
totleng ~ dnorm(mu,sigma),
mu <- alpha +
alpha_acs[acs_id] +
beta_sp[species_id] +
beta_shade*shade +
beta_sp_shade[species_id]*shade,
alpha ~ dnorm(53,20),
alpha_acs[acs_id] ~ dnorm(0,sigma_acs),
beta_sp[species_id] ~ dnorm(0,10),
beta_shade ~ dnorm(0,20),
beta_sp_shade[species_id] ~ dnorm(0,10),
c(sigma,sigma_acs) ~ dcauchy(0,1)
),
data = tomato.small,
iter = 4000, # Rhat a bit high with defaults
warmup = 1000,
chains=4,
cores=2
)##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 4e-06 seconds (Warm-up)
## 0.000458 seconds (Sampling)
## 0.000462 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 1200 / 12000 ]
[ 2400 / 12000 ]
[ 3600 / 12000 ]
[ 4800 / 12000 ]
[ 6000 / 12000 ]
[ 7200 / 12000 ]
[ 8400 / 12000 ]
[ 9600 / 12000 ]
[ 10800 / 12000 ]
[ 12000 / 12000 ]plot(m.shade.species.int.acs,ask=FALSE)## Waiting to draw page 2 of 4
## Waiting to draw page 3 of 4
## Waiting to draw page 4 of 4
par(mfrow=c(1,1))
precis(m.shade.species.int.acs,depth = 2)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha 41.51 4.45 34.40 48.48 2958 1
## alpha_acs[1] -5.29 3.67 -10.92 0.71 5211 1
## alpha_acs[2] -10.62 3.51 -16.10 -4.89 4596 1
## alpha_acs[3] -3.46 3.37 -8.82 1.92 4256 1
## alpha_acs[4] 3.25 3.39 -2.42 8.45 4536 1
## alpha_acs[5] 13.65 3.43 8.33 19.23 4113 1
## alpha_acs[6] 2.59 3.40 -2.82 7.97 4251 1
## alpha_acs[7] 9.84 3.47 4.46 15.48 4689 1
## alpha_acs[8] 6.43 3.28 0.95 11.46 4880 1
## alpha_acs[9] 3.54 3.31 -2.03 8.58 4603 1
## alpha_acs[10] 3.92 3.91 -2.37 10.11 6203 1
## alpha_acs[11] 2.26 3.37 -3.40 7.33 4482 1
## alpha_acs[12] -9.69 3.45 -15.03 -4.01 4875 1
## alpha_acs[13] -3.96 3.65 -9.74 1.80 5408 1
## alpha_acs[14] -0.62 3.53 -6.39 4.92 5236 1
## alpha_acs[15] -3.49 4.92 -10.88 4.86 8056 1
## alpha_acs[16] 1.30 3.29 -4.28 6.35 4206 1
## alpha_acs[17] 4.76 3.36 -0.59 10.18 4948 1
## alpha_acs[18] -12.05 3.45 -17.56 -6.58 5163 1
## alpha_acs[19] 14.99 3.49 9.38 20.45 4413 1
## alpha_acs[20] -6.10 3.23 -11.34 -1.09 4757 1
## alpha_acs[21] 1.09 3.25 -4.15 6.27 4622 1
## alpha_acs[22] 6.27 4.45 -0.61 13.48 6648 1
## alpha_acs[23] -2.14 3.77 -8.11 3.89 5467 1
## alpha_acs[24] -4.20 3.33 -9.47 1.11 4277 1
## alpha_acs[25] -5.02 3.30 -10.13 0.36 4121 1
## alpha_acs[26] -0.30 3.58 -6.03 5.39 5309 1
## alpha_acs[27] 1.84 3.29 -3.55 6.98 4846 1
## alpha_acs[28] -4.60 3.42 -10.23 0.66 5107 1
## alpha_acs[29] 1.17 3.46 -4.30 6.78 5227 1
## alpha_acs[30] 6.80 3.34 1.50 12.10 4767 1
## alpha_acs[31] 10.41 3.42 4.97 15.85 5080 1
## alpha_acs[32] 0.58 3.62 -5.08 6.47 5236 1
## alpha_acs[33] -3.24 3.59 -8.88 2.53 5262 1
## alpha_acs[34] -8.66 3.42 -13.98 -3.12 4760 1
## alpha_acs[35] -5.58 3.60 -11.34 0.09 5101 1
## alpha_acs[36] -7.12 4.38 -13.77 0.09 6366 1
## beta_sp[1] -4.02 4.96 -12.02 3.80 3036 1
## beta_sp[2] 2.45 4.93 -5.53 10.20 2927 1
## beta_sp[3] 2.57 5.05 -5.53 10.50 3228 1
## beta_sp[4] -10.77 5.03 -18.55 -2.75 3424 1
## beta_sp[5] 6.26 4.98 -1.65 14.29 2884 1
## beta_shade 22.07 4.32 15.05 28.81 2714 1
## beta_sp_shade[1] 5.70 4.54 -1.41 12.99 2908 1
## beta_sp_shade[2] -5.14 4.57 -12.57 2.08 2975 1
## beta_sp_shade[3] -1.91 4.54 -8.86 5.63 2985 1
## beta_sp_shade[4] 5.08 4.65 -2.39 12.37 3148 1
## beta_sp_shade[5] 2.37 4.55 -4.92 9.76 2949 1
## sigma 13.34 0.30 12.84 13.81 8036 1
## sigma_acs 7.43 1.07 5.73 9.05 6140 1compare(m.shade.species.int,m.shade.species.int.acs)## WAIC pWAIC dWAIC weight SE dSE
## m.shade.species.int.acs 8123.1 38.1 0.0 1 55.18 NA
## m.shade.species.int 8335.2 11.0 212.1 0 53.38 28.15coeftab(m.shade.species.int,m.shade.species.int.acs)## m.shade.species.int m.shade.species.int.acs
## alpha 41.62 41.51
## beta_sp[1] -4.31 -4.02
## beta_sp[2] 2.62 2.45
## beta_sp[3] 3.10 2.57
## beta_sp[4] -11.82 -10.77
## beta_sp[5] 7.05 6.26
## beta_shade 22.56 22.07
## beta_sp_shade[1] 5.51 5.70
## beta_sp_shade[2] -5.41 -5.14
## beta_sp_shade[3] -2.04 -1.91
## beta_sp_shade[4] 4.01 5.08
## beta_sp_shade[5] 2.06 2.37
## sigma 15.03 13.34
## alpha_acs[1] NA -5.29
## alpha_acs[2] NA -10.62
## alpha_acs[3] NA -3.46
## alpha_acs[4] NA 3.25
## alpha_acs[5] NA 13.65
## alpha_acs[6] NA 2.59
## alpha_acs[7] NA 9.84
## alpha_acs[8] NA 6.43
## alpha_acs[9] NA 3.54
## alpha_acs[10] NA 3.92
## alpha_acs[11] NA 2.26
## alpha_acs[12] NA -9.69
## alpha_acs[13] NA -3.96
## alpha_acs[14] NA -0.62
## alpha_acs[15] NA -3.49
## alpha_acs[16] NA 1.3
## alpha_acs[17] NA 4.76
## alpha_acs[18] NA -12.05
## alpha_acs[19] NA 14.99
## alpha_acs[20] NA -6.1
## alpha_acs[21] NA 1.09
## alpha_acs[22] NA 6.27
## alpha_acs[23] NA -2.14
## alpha_acs[24] NA -4.2
## alpha_acs[25] NA -5.02
## alpha_acs[26] NA -0.3
## alpha_acs[27] NA 1.84
## alpha_acs[28] NA -4.6
## alpha_acs[29] NA 1.17
## alpha_acs[30] NA 6.8
## alpha_acs[31] NA 10.41
## alpha_acs[32] NA 0.58
## alpha_acs[33] NA -3.24
## alpha_acs[34] NA -8.66
## alpha_acs[35] NA -5.58
## alpha_acs[36] NA -7.12
## sigma_acs NA 7.43
## nobs 1008 1008plot(coeftab(m.shade.species.int,m.shade.species.int.acs),pars=1:13)
Adding acs does not change most of the shared coefficients by much, although it drops sigma.
Add shelf random effect
table(tomato.small$shade, tomato.small$shelf_id)##
## 1 2 3 4 5 6
## 0 0 0 0 174 125 196
## 1 161 174 178 0 0 0table(tomato.small$acs_id,tomato.small$shelf_id)##
## 1 2 3 4 5 6
## 1 0 5 5 0 4 6
## 2 2 6 6 3 2 6
## 3 6 6 6 6 4 6
## 4 6 6 6 6 3 6
## 5 6 6 6 5 5 6
## 6 5 5 5 5 5 6
## 7 5 6 6 6 2 6
## 8 7 4 6 9 5 6
## 9 6 5 6 6 5 6
## 10 0 6 2 0 3 5
## 11 6 6 6 9 4 5
## 12 2 4 6 3 5 6
## 13 3 6 5 1 3 6
## 14 2 2 6 3 4 6
## 15 0 2 0 0 0 3
## 16 8 6 6 9 5 6
## 17 6 5 5 7 4 6
## 18 4 6 3 5 4 6
## 19 5 6 5 5 3 6
## 20 10 2 6 11 3 5
## 21 8 4 4 9 5 6
## 22 0 5 0 0 0 3
## 23 3 4 5 0 3 6
## 24 6 6 6 6 5 6
## 25 8 6 6 9 4 6
## 26 3 3 6 4 2 4
## 27 7 6 5 9 4 6
## 28 4 6 5 5 4 5
## 29 4 4 5 5 3 5
## 30 5 6 7 5 3 6
## 31 6 3 6 4 4 5
## 32 3 4 4 6 3 5
## 33 5 5 6 5 1 6
## 34 5 6 5 5 4 6
## 35 3 5 3 3 4 6
## 36 2 1 3 0 3 1most accessions are on most shelves
m.shade.species.int.acs.shelf <- map2stan(
alist(
totleng ~ dnorm(mu,sigma),
mu <- alpha +
alpha_acs[acs_id] +
alpha_shelf[shelf_id] +
beta_sp[species_id] +
beta_shade*shade +
beta_sp_shade[species_id]*shade,
alpha ~ dnorm(53,20),
alpha_acs[acs_id] ~ dnorm(0,sigma_acs),
alpha_shelf[shelf_id] ~ dnorm(0,sigma_shelf),
beta_sp[species_id] ~ dnorm(0,10),
beta_shade ~ dnorm(0,20),
beta_sp_shade[species_id] ~ dnorm(0,10),
c(sigma,sigma_acs) ~ dcauchy(0,1),
sigma_shelf ~ dexp(1) # with only 3 shelves per trt might be hard to estimate w. cauchy
),
data = tomato.small,
iter = 4000, # Rhat a bit high with defaults
warmup = 1000,
chains=4,
cores=2
)##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 5e-06 seconds (Warm-up)
## 0.000539 seconds (Sampling)
## 0.000544 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 1200 / 12000 ]
[ 2400 / 12000 ]
[ 3600 / 12000 ]
[ 4800 / 12000 ]
[ 6000 / 12000 ]
[ 7200 / 12000 ]
[ 8400 / 12000 ]
[ 9600 / 12000 ]
[ 10800 / 12000 ]
[ 12000 / 12000 ]save.image(file="tomato1212.Rdata")plot(m.shade.species.int.acs.shelf,ask=FALSE)## Waiting to draw page 2 of 4
## Waiting to draw page 3 of 4
## Waiting to draw page 4 of 4
par(mfrow=c(1,1))
precis(m.shade.species.int.acs.shelf,depth = 2)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha 41.19 4.64 33.47 48.30 4939 1
## alpha_acs[1] -6.46 3.71 -12.54 -0.61 7764 1
## alpha_acs[2] -11.28 3.50 -16.92 -5.75 6628 1
## alpha_acs[3] -3.36 3.40 -8.67 2.17 6679 1
## alpha_acs[4] 3.41 3.40 -2.08 8.83 6843 1
## alpha_acs[5] 13.85 3.41 8.40 19.37 6813 1
## alpha_acs[6] 2.73 3.40 -2.84 8.05 6354 1
## alpha_acs[7] 9.92 3.39 4.44 15.19 6943 1
## alpha_acs[8] 6.94 3.31 1.57 12.17 6281 1
## alpha_acs[9] 3.69 3.35 -1.60 8.92 6709 1
## alpha_acs[10] 2.99 3.89 -3.29 9.08 8034 1
## alpha_acs[11] 2.71 3.29 -2.44 8.06 6086 1
## alpha_acs[12] -10.11 3.49 -15.63 -4.56 6961 1
## alpha_acs[13] -4.08 3.61 -9.61 1.91 6835 1
## alpha_acs[14] -0.78 3.55 -6.33 4.92 7230 1
## alpha_acs[15] -5.43 4.90 -13.16 2.45 9651 1
## alpha_acs[16] 1.86 3.26 -3.42 6.94 6496 1
## alpha_acs[17] 5.13 3.49 -0.16 10.97 5678 1
## alpha_acs[18] -12.23 3.56 -18.02 -6.66 5534 1
## alpha_acs[19] 15.09 3.44 9.55 20.40 6176 1
## alpha_acs[20] -5.10 3.29 -10.36 0.14 6053 1
## alpha_acs[21] 1.64 3.32 -3.46 7.09 6287 1
## alpha_acs[22] 4.55 4.46 -2.40 11.78 8149 1
## alpha_acs[23] -2.15 3.68 -7.81 3.99 7161 1
## alpha_acs[24] -4.07 3.37 -9.68 1.08 6767 1
## alpha_acs[25] -4.57 3.35 -9.72 0.90 6407 1
## alpha_acs[26] -0.08 3.68 -5.87 5.82 5986 1
## alpha_acs[27] 2.30 3.43 -3.26 7.70 5488 1
## alpha_acs[28] -4.60 3.53 -10.12 1.01 5919 1
## alpha_acs[29] 1.37 3.57 -4.57 6.82 6222 1
## alpha_acs[30] 6.84 3.46 1.15 12.19 5446 1
## alpha_acs[31] 10.79 3.45 5.31 16.23 6680 1
## alpha_acs[32] 1.29 3.59 -4.41 6.96 7396 1
## alpha_acs[33] -2.71 3.57 -8.52 2.85 7366 1
## alpha_acs[34] -8.68 3.34 -14.17 -3.64 7032 1
## alpha_acs[35] -5.99 3.50 -11.74 -0.53 7064 1
## alpha_acs[36] -6.42 4.37 -13.42 0.58 7852 1
## alpha_shelf[1] -3.09 1.87 -6.02 -0.21 6722 1
## alpha_shelf[2] 3.62 1.86 0.62 6.43 5914 1
## alpha_shelf[3] -0.15 1.84 -3.15 2.61 6508 1
## alpha_shelf[4] -2.37 1.88 -5.32 0.55 6682 1
## alpha_shelf[5] -1.05 1.87 -4.05 1.78 7218 1
## alpha_shelf[6] 2.70 1.85 -0.20 5.59 7240 1
## beta_sp[1] -3.72 4.94 -12.02 3.80 5048 1
## beta_sp[2] 2.74 4.97 -5.17 10.79 5171 1
## beta_sp[3] 2.81 4.96 -5.29 10.53 4961 1
## beta_sp[4] -11.21 5.04 -19.54 -3.50 5380 1
## beta_sp[5] 6.58 4.96 -1.21 14.67 5274 1
## beta_shade 22.00 4.91 14.18 29.88 4781 1
## beta_sp_shade[1] 5.40 4.61 -1.91 12.84 5026 1
## beta_sp_shade[2] -5.33 4.58 -12.71 1.85 4661 1
## beta_sp_shade[3] -1.69 4.58 -8.94 5.60 5148 1
## beta_sp_shade[4] 5.05 4.70 -2.72 12.36 4875 1
## beta_sp_shade[5] 2.03 4.60 -5.44 9.21 5084 1
## sigma 13.05 0.30 12.58 13.53 9461 1
## sigma_acs 7.54 1.09 5.90 9.25 6836 1
## sigma_shelf 2.85 0.89 1.53 4.08 6032 1compare(m.shade,m.shade.species,m.shade.species.int,m.shade.species.int.acs,m.shade.species.int.acs.shelf)## WAIC pWAIC dWAIC weight SE dSE
## m.shade.species.int.acs.shelf 8083.3 42.2 0.0 1 56.41 NA
## m.shade.species.int.acs 8123.1 38.1 39.8 0 55.18 12.03
## m.shade.species.int 8335.2 11.0 251.9 0 53.38 30.17
## m.shade.species 8348.2 6.9 264.9 0 53.08 31.11
## m.shade 8458.3 3.1 374.9 0 52.60 38.43coeftab(m.shade.species.int.acs,m.shade.species.int.acs.shelf)## m.shade.species.int.acs m.shade.species.int.acs.shelf
## alpha 41.51 41.19
## alpha_acs[1] -5.29 -6.46
## alpha_acs[2] -10.62 -11.28
## alpha_acs[3] -3.46 -3.36
## alpha_acs[4] 3.25 3.41
## alpha_acs[5] 13.65 13.85
## alpha_acs[6] 2.59 2.73
## alpha_acs[7] 9.84 9.92
## alpha_acs[8] 6.43 6.94
## alpha_acs[9] 3.54 3.69
## alpha_acs[10] 3.92 2.99
## alpha_acs[11] 2.26 2.71
## alpha_acs[12] -9.69 -10.11
## alpha_acs[13] -3.96 -4.08
## alpha_acs[14] -0.62 -0.78
## alpha_acs[15] -3.49 -5.43
## alpha_acs[16] 1.30 1.86
## alpha_acs[17] 4.76 5.13
## alpha_acs[18] -12.05 -12.23
## alpha_acs[19] 14.99 15.09
## alpha_acs[20] -6.1 -5.1
## alpha_acs[21] 1.09 1.64
## alpha_acs[22] 6.27 4.55
## alpha_acs[23] -2.14 -2.15
## alpha_acs[24] -4.20 -4.07
## alpha_acs[25] -5.02 -4.57
## alpha_acs[26] -0.30 -0.08
## alpha_acs[27] 1.84 2.30
## alpha_acs[28] -4.6 -4.6
## alpha_acs[29] 1.17 1.37
## alpha_acs[30] 6.80 6.84
## alpha_acs[31] 10.41 10.79
## alpha_acs[32] 0.58 1.29
## alpha_acs[33] -3.24 -2.71
## alpha_acs[34] -8.66 -8.68
## alpha_acs[35] -5.58 -5.99
## alpha_acs[36] -7.12 -6.42
## beta_sp[1] -4.02 -3.72
## beta_sp[2] 2.45 2.74
## beta_sp[3] 2.57 2.81
## beta_sp[4] -10.77 -11.21
## beta_sp[5] 6.26 6.58
## beta_shade 22.07 22.00
## beta_sp_shade[1] 5.7 5.4
## beta_sp_shade[2] -5.14 -5.33
## beta_sp_shade[3] -1.91 -1.69
## beta_sp_shade[4] 5.08 5.05
## beta_sp_shade[5] 2.37 2.03
## sigma 13.34 13.05
## sigma_acs 7.43 7.54
## alpha_shelf[1] NA -3.09
## alpha_shelf[2] NA 3.62
## alpha_shelf[3] NA -0.15
## alpha_shelf[4] NA -2.37
## alpha_shelf[5] NA -1.05
## alpha_shelf[6] NA 2.7
## sigma_shelf NA 2.85
## nobs 1008 1008plot(coeftab(m.shade.species.int.acs,m.shade.species.int.acs.shelf),pars=c(1,38:50))
including shelf does not change the estimates much, although that model is apparently preferred
reparamterize with no “main” shade effect
m.shade.species.sep.acs.shelf <- map2stan(
alist(
totleng ~ dnorm(mu,sigma),
mu <- alpha +
alpha_acs[acs_id] +
alpha_shelf[shelf_id] +
beta_sp[species_id] +
beta_sp_shade[species_id]*shade,
alpha ~ dnorm(53,20),
alpha_acs[acs_id] ~ dnorm(0,sigma_acs),
alpha_shelf[shelf_id] ~ dnorm(0,sigma_shelf),
beta_sp[species_id] ~ dnorm(0,10),
beta_sp_shade[species_id] ~ dnorm(0,20),
c(sigma,sigma_acs) ~ dcauchy(0,1),
sigma_shelf ~ dexp(1) # with only 3 shelves per trt might be hard to estimate w. cauchy
),
data = tomato.small,
iter = 4000, # Rhat a bit high with defaults
warmup = 1000,
chains=4,
cores=4
)##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 3e-06 seconds (Warm-up)
## 0.000299 seconds (Sampling)
## 0.000302 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 1200 / 12000 ]
[ 2400 / 12000 ]
[ 3600 / 12000 ]
[ 4800 / 12000 ]
[ 6000 / 12000 ]
[ 7200 / 12000 ]
[ 8400 / 12000 ]
[ 9600 / 12000 ]
[ 10800 / 12000 ]
[ 12000 / 12000 ]save.image(file="tomato1212.Rdata")plot(m.shade.species.sep.acs.shelf,ask=FALSE)## Waiting to draw page 2 of 4
## Waiting to draw page 3 of 4
## Waiting to draw page 4 of 4
par(mfrow=c(1,1))
precis(m.shade.species.sep.acs.shelf,depth = 2)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha 41.82 4.84 34.10 49.49 2074 1
## alpha_acs[1] -6.55 3.68 -12.17 -0.52 3673 1
## alpha_acs[2] -11.24 3.53 -17.00 -5.77 3665 1
## alpha_acs[3] -3.37 3.38 -8.51 2.21 2942 1
## alpha_acs[4] 3.34 3.40 -2.20 8.69 3031 1
## alpha_acs[5] 13.83 3.44 8.42 19.35 3121 1
## alpha_acs[6] 2.66 3.46 -3.06 7.97 3343 1
## alpha_acs[7] 10.04 3.46 4.66 15.64 3337 1
## alpha_acs[8] 6.90 3.27 1.87 12.30 3039 1
## alpha_acs[9] 3.70 3.31 -1.49 9.04 3093 1
## alpha_acs[10] 3.01 3.87 -3.26 9.08 4292 1
## alpha_acs[11] 2.78 3.34 -2.47 8.05 3230 1
## alpha_acs[12] -10.13 3.49 -15.69 -4.60 3200 1
## alpha_acs[13] -4.02 3.70 -10.05 1.75 3979 1
## alpha_acs[14] -0.78 3.49 -6.43 4.65 3570 1
## alpha_acs[15] -5.42 4.87 -13.24 2.31 6145 1
## alpha_acs[16] 1.94 3.30 -3.12 7.33 3237 1
## alpha_acs[17] 5.08 3.43 -0.41 10.52 2773 1
## alpha_acs[18] -12.32 3.51 -17.97 -6.69 2696 1
## alpha_acs[19] 15.14 3.48 9.43 20.42 3409 1
## alpha_acs[20] -5.08 3.25 -10.10 0.25 2945 1
## alpha_acs[21] 1.63 3.28 -3.78 6.56 3090 1
## alpha_acs[22] 4.63 4.47 -2.06 12.29 5431 1
## alpha_acs[23] -2.12 3.75 -8.07 3.90 4249 1
## alpha_acs[24] -4.08 3.39 -9.44 1.41 3111 1
## alpha_acs[25] -4.58 3.32 -9.83 0.65 2954 1
## alpha_acs[26] -0.17 3.68 -5.99 5.69 2989 1
## alpha_acs[27] 2.22 3.36 -2.88 7.74 2909 1
## alpha_acs[28] -4.70 3.49 -10.27 0.93 2659 1
## alpha_acs[29] 1.30 3.53 -4.28 6.96 2887 1
## alpha_acs[30] 6.79 3.43 1.50 12.43 2745 1
## alpha_acs[31] 10.80 3.39 5.44 16.26 3152 1
## alpha_acs[32] 1.34 3.64 -4.54 7.06 3974 1
## alpha_acs[33] -2.64 3.60 -8.03 3.45 3908 1
## alpha_acs[34] -8.60 3.42 -14.13 -3.26 3260 1
## alpha_acs[35] -5.92 3.53 -11.66 -0.41 3504 1
## alpha_acs[36] -6.31 4.39 -13.16 0.75 5423 1
## alpha_shelf[1] -2.40 1.95 -5.54 0.52 1998 1
## alpha_shelf[2] 4.39 2.06 1.23 7.54 1762 1
## alpha_shelf[3] 0.57 1.96 -2.51 3.52 1845 1
## alpha_shelf[4] -3.09 2.02 -6.31 -0.02 2442 1
## alpha_shelf[5] -1.75 2.02 -4.88 1.41 2557 1
## alpha_shelf[6] 2.02 1.94 -0.96 5.07 2599 1
## beta_sp[1] -3.57 5.11 -11.67 4.70 2196 1
## beta_sp[2] 2.98 5.08 -4.95 11.10 2163 1
## beta_sp[3] 2.98 5.09 -5.43 10.83 2256 1
## beta_sp[4] -11.17 5.21 -19.83 -3.13 2243 1
## beta_sp[5] 6.62 5.08 -1.59 14.69 2051 1
## beta_sp_shade[1] 25.96 3.18 21.28 31.18 1702 1
## beta_sp_shade[2] 14.97 3.13 10.27 20.07 1740 1
## beta_sp_shade[3] 18.67 3.13 13.98 23.66 1738 1
## beta_sp_shade[4] 25.55 3.44 20.06 30.81 1903 1
## beta_sp_shade[5] 22.51 3.16 17.64 27.55 1765 1
## sigma 13.04 0.29 12.55 13.49 8658 1
## sigma_acs 7.53 1.08 5.81 9.15 4816 1
## sigma_shelf 2.98 0.96 1.51 4.29 2904 1compare(m.shade,m.shade.species,m.shade.species.int,m.shade.species.int.acs,m.shade.species.int.acs.shelf,m.shade.species.sep.acs.shelf)## WAIC pWAIC dWAIC weight SE dSE
## m.shade.species.int.acs.shelf 8083.3 42.2 0.0 0.54 56.41 NA
## m.shade.species.sep.acs.shelf 8083.6 42.3 0.3 0.46 56.44 0.44
## m.shade.species.int.acs 8123.1 38.1 39.8 0.00 55.18 12.03
## m.shade.species.int 8335.2 11.0 251.9 0.00 53.38 30.17
## m.shade.species 8348.2 6.9 264.9 0.00 53.08 31.11
## m.shade 8458.3 3.1 374.9 0.00 52.60 38.43coeftab(m.shade.species.sep.acs.shelf,m.shade.species.int.acs.shelf)## m.shade.species.sep.acs.shelf
## alpha 41.82
## alpha_acs[1] -6.55
## alpha_acs[2] -11.24
## alpha_acs[3] -3.37
## alpha_acs[4] 3.34
## alpha_acs[5] 13.83
## alpha_acs[6] 2.66
## alpha_acs[7] 10.04
## alpha_acs[8] 6.90
## alpha_acs[9] 3.70
## alpha_acs[10] 3.01
## alpha_acs[11] 2.78
## alpha_acs[12] -10.13
## alpha_acs[13] -4.02
## alpha_acs[14] -0.78
## alpha_acs[15] -5.42
## alpha_acs[16] 1.94
## alpha_acs[17] 5.08
## alpha_acs[18] -12.32
## alpha_acs[19] 15.14
## alpha_acs[20] -5.08
## alpha_acs[21] 1.63
## alpha_acs[22] 4.63
## alpha_acs[23] -2.12
## alpha_acs[24] -4.08
## alpha_acs[25] -4.58
## alpha_acs[26] -0.17
## alpha_acs[27] 2.22
## alpha_acs[28] -4.7
## alpha_acs[29] 1.30
## alpha_acs[30] 6.79
## alpha_acs[31] 10.80
## alpha_acs[32] 1.34
## alpha_acs[33] -2.64
## alpha_acs[34] -8.60
## alpha_acs[35] -5.92
## alpha_acs[36] -6.31
## alpha_shelf[1] -2.40
## alpha_shelf[2] 4.39
## alpha_shelf[3] 0.57
## alpha_shelf[4] -3.09
## alpha_shelf[5] -1.75
## alpha_shelf[6] 2.02
## beta_sp[1] -3.57
## beta_sp[2] 2.98
## beta_sp[3] 2.98
## beta_sp[4] -11.17
## beta_sp[5] 6.62
## beta_sp_shade[1] 25.96
## beta_sp_shade[2] 14.97
## beta_sp_shade[3] 18.67
## beta_sp_shade[4] 25.55
## beta_sp_shade[5] 22.51
## sigma 13.04
## sigma_acs 7.53
## sigma_shelf 2.98
## beta_shade NA
## nobs 1008
## m.shade.species.int.acs.shelf
## alpha 41.19
## alpha_acs[1] -6.46
## alpha_acs[2] -11.28
## alpha_acs[3] -3.36
## alpha_acs[4] 3.41
## alpha_acs[5] 13.85
## alpha_acs[6] 2.73
## alpha_acs[7] 9.92
## alpha_acs[8] 6.94
## alpha_acs[9] 3.69
## alpha_acs[10] 2.99
## alpha_acs[11] 2.71
## alpha_acs[12] -10.11
## alpha_acs[13] -4.08
## alpha_acs[14] -0.78
## alpha_acs[15] -5.43
## alpha_acs[16] 1.86
## alpha_acs[17] 5.13
## alpha_acs[18] -12.23
## alpha_acs[19] 15.09
## alpha_acs[20] -5.10
## alpha_acs[21] 1.64
## alpha_acs[22] 4.55
## alpha_acs[23] -2.15
## alpha_acs[24] -4.07
## alpha_acs[25] -4.57
## alpha_acs[26] -0.08
## alpha_acs[27] 2.30
## alpha_acs[28] -4.6
## alpha_acs[29] 1.37
## alpha_acs[30] 6.84
## alpha_acs[31] 10.79
## alpha_acs[32] 1.29
## alpha_acs[33] -2.71
## alpha_acs[34] -8.68
## alpha_acs[35] -5.99
## alpha_acs[36] -6.42
## alpha_shelf[1] -3.09
## alpha_shelf[2] 3.62
## alpha_shelf[3] -0.15
## alpha_shelf[4] -2.37
## alpha_shelf[5] -1.05
## alpha_shelf[6] 2.70
## beta_sp[1] -3.72
## beta_sp[2] 2.74
## beta_sp[3] 2.81
## beta_sp[4] -11.21
## beta_sp[5] 6.58
## beta_sp_shade[1] 5.40
## beta_sp_shade[2] -5.33
## beta_sp_shade[3] -1.69
## beta_sp_shade[4] 5.05
## beta_sp_shade[5] 2.03
## sigma 13.05
## sigma_acs 7.54
## sigma_shelf 2.85
## beta_shade 22
## nobs 1008plot(coeftab(m.shade.species.sep.acs.shelf))
This parameterization leads to the same WAIC, so it is an equivalent model. Interestingly though, the SD on the shade coefficients is smaller, presumably because one fewer parameter is being estimated.
Can we do no main alpha?
m.shade.species.sep.acs.shelf.no.alpha <- map2stan(
alist(
totleng ~ dnorm(mu,sigma),
mu <-
alpha_acs[acs_id] +
alpha_shelf[shelf_id] +
beta_sp[species_id] +
beta_sp_shade[species_id]*shade,
alpha_acs[acs_id] ~ dnorm(0,sigma_acs),
alpha_shelf[shelf_id] ~ dnorm(0,sigma_shelf),
beta_sp[species_id] ~ dnorm(42,20),
beta_sp_shade[species_id] ~ dnorm(0,20),
c(sigma,sigma_acs) ~ dcauchy(0,1),
sigma_shelf ~ dexp(1) # with only 3 shelves per trt might be hard to estimate w. cauchy
),
data = tomato.small,
iter = 4000, # Rhat a bit high with defaults
warmup = 1000,
chains=4,
cores=4
)##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 5e-06 seconds (Warm-up)
## 0.000256 seconds (Sampling)
## 0.000261 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 1200 / 12000 ]
[ 2400 / 12000 ]
[ 3600 / 12000 ]
[ 4800 / 12000 ]
[ 6000 / 12000 ]
[ 7200 / 12000 ]
[ 8400 / 12000 ]
[ 9600 / 12000 ]
[ 10800 / 12000 ]
[ 12000 / 12000 ]save.image(file="tomato1212.Rdata")plot(m.shade.species.sep.acs.shelf.no.alpha,ask=FALSE)## Waiting to draw page 2 of 4
## Waiting to draw page 3 of 4
## Waiting to draw page 4 of 4
par(mfrow=c(1,1))
precis(m.shade.species.sep.acs.shelf.no.alpha,depth = 2)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha_acs[1] -6.56 3.79 -12.80 -0.66 3211 1
## alpha_acs[2] -11.70 3.58 -17.41 -5.99 3545 1
## alpha_acs[3] -3.38 3.48 -9.11 2.03 2588 1
## alpha_acs[4] 3.34 3.51 -2.25 8.91 2646 1
## alpha_acs[5] 13.82 3.52 8.36 19.60 2530 1
## alpha_acs[6] 2.64 3.52 -2.99 8.12 2577 1
## alpha_acs[7] 9.60 3.48 3.98 15.08 3175 1
## alpha_acs[8] 6.77 3.34 1.44 12.04 2967 1
## alpha_acs[9] 3.54 3.38 -1.77 9.05 3025 1
## alpha_acs[10] 3.60 4.06 -2.95 10.04 3570 1
## alpha_acs[11] 2.32 3.39 -3.24 7.61 3084 1
## alpha_acs[12] -10.32 3.53 -15.84 -4.61 3312 1
## alpha_acs[13] -3.39 3.77 -9.49 2.37 3266 1
## alpha_acs[14] -0.94 3.57 -6.52 4.81 3506 1
## alpha_acs[15] -5.62 5.03 -14.00 2.13 6728 1
## alpha_acs[16] 1.46 3.37 -4.01 6.68 2963 1
## alpha_acs[17] 5.19 3.51 -0.41 10.83 3115 1
## alpha_acs[18] -12.21 3.61 -18.08 -6.56 3180 1
## alpha_acs[19] 14.72 3.50 9.14 20.28 3231 1
## alpha_acs[20] -5.28 3.33 -10.54 0.08 2931 1
## alpha_acs[21] 1.46 3.34 -3.79 6.87 3072 1
## alpha_acs[22] 5.18 4.58 -2.27 12.20 5099 1
## alpha_acs[23] -1.51 3.83 -7.88 4.27 3224 1
## alpha_acs[24] -4.11 3.51 -9.44 1.69 2543 1
## alpha_acs[25] -4.64 3.44 -10.42 0.61 2500 1
## alpha_acs[26] -0.02 3.76 -6.25 5.81 3440 1
## alpha_acs[27] 2.34 3.48 -3.19 7.92 3042 1
## alpha_acs[28] -4.56 3.61 -10.43 1.02 3173 1
## alpha_acs[29] 1.42 3.62 -4.33 7.16 3275 1
## alpha_acs[30] 6.90 3.52 1.32 12.60 3098 1
## alpha_acs[31] 10.64 3.51 4.82 16.06 3210 1
## alpha_acs[32] 1.98 3.76 -3.92 8.11 3307 1
## alpha_acs[33] -1.99 3.67 -7.92 3.82 3229 1
## alpha_acs[34] -9.05 3.44 -14.80 -3.84 3069 1
## alpha_acs[35] -6.39 3.63 -12.31 -0.76 3341 1
## alpha_acs[36] -5.79 4.39 -12.96 1.02 4407 1
## alpha_shelf[1] -2.40 1.93 -5.60 0.48 3024 1
## alpha_shelf[2] 4.37 2.04 1.18 7.48 2615 1
## alpha_shelf[3] 0.57 1.95 -2.48 3.54 2686 1
## alpha_shelf[4] -3.03 1.95 -6.01 0.12 1907 1
## alpha_shelf[5] -1.68 1.97 -4.79 1.38 2037 1
## alpha_shelf[6] 2.11 1.87 -0.72 5.18 2002 1
## beta_sp[1] 38.01 3.58 32.31 43.74 2052 1
## beta_sp[2] 44.79 3.49 39.30 50.52 1742 1
## beta_sp[3] 44.91 3.42 38.93 49.92 1943 1
## beta_sp[4] 29.75 3.78 23.62 35.62 2173 1
## beta_sp[5] 48.92 3.49 43.31 54.37 2141 1
## beta_sp_shade[1] 26.07 3.17 21.35 31.21 1957 1
## beta_sp_shade[2] 15.01 3.11 10.19 19.93 1869 1
## beta_sp_shade[3] 18.76 3.14 13.94 23.68 1971 1
## beta_sp_shade[4] 25.85 3.49 20.34 31.34 2245 1
## beta_sp_shade[5] 22.48 3.16 17.59 27.42 1966 1
## sigma 13.04 0.30 12.57 13.52 11087 1
## sigma_acs 7.59 1.08 5.96 9.33 5488 1
## sigma_shelf 2.97 0.95 1.59 4.32 3635 1compare(m.shade,m.shade.species,m.shade.species.int,m.shade.species.int.acs,m.shade.species.int.acs.shelf,m.shade.species.sep.acs.shelf,m.shade.species.sep.acs.shelf.no.alpha)## WAIC pWAIC dWAIC weight SE
## m.shade.species.int.acs.shelf 8083.3 42.2 0.0 0.36 56.41
## m.shade.species.sep.acs.shelf.no.alpha 8083.5 42.4 0.2 0.32 56.47
## m.shade.species.sep.acs.shelf 8083.6 42.3 0.3 0.31 56.44
## m.shade.species.int.acs 8123.1 38.1 39.8 0.00 55.18
## m.shade.species.int 8335.2 11.0 251.9 0.00 53.38
## m.shade.species 8348.2 6.9 264.9 0.00 53.08
## m.shade 8458.3 3.1 374.9 0.00 52.60
## dSE
## m.shade.species.int.acs.shelf NA
## m.shade.species.sep.acs.shelf.no.alpha 0.50
## m.shade.species.sep.acs.shelf 0.44
## m.shade.species.int.acs 12.03
## m.shade.species.int 30.17
## m.shade.species 31.11
## m.shade 38.43coeftab(m.shade.species.sep.acs.shelf,m.shade.species.int.acs.shelf)## m.shade.species.sep.acs.shelf
## alpha 41.82
## alpha_acs[1] -6.55
## alpha_acs[2] -11.24
## alpha_acs[3] -3.37
## alpha_acs[4] 3.34
## alpha_acs[5] 13.83
## alpha_acs[6] 2.66
## alpha_acs[7] 10.04
## alpha_acs[8] 6.90
## alpha_acs[9] 3.70
## alpha_acs[10] 3.01
## alpha_acs[11] 2.78
## alpha_acs[12] -10.13
## alpha_acs[13] -4.02
## alpha_acs[14] -0.78
## alpha_acs[15] -5.42
## alpha_acs[16] 1.94
## alpha_acs[17] 5.08
## alpha_acs[18] -12.32
## alpha_acs[19] 15.14
## alpha_acs[20] -5.08
## alpha_acs[21] 1.63
## alpha_acs[22] 4.63
## alpha_acs[23] -2.12
## alpha_acs[24] -4.08
## alpha_acs[25] -4.58
## alpha_acs[26] -0.17
## alpha_acs[27] 2.22
## alpha_acs[28] -4.7
## alpha_acs[29] 1.30
## alpha_acs[30] 6.79
## alpha_acs[31] 10.80
## alpha_acs[32] 1.34
## alpha_acs[33] -2.64
## alpha_acs[34] -8.60
## alpha_acs[35] -5.92
## alpha_acs[36] -6.31
## alpha_shelf[1] -2.40
## alpha_shelf[2] 4.39
## alpha_shelf[3] 0.57
## alpha_shelf[4] -3.09
## alpha_shelf[5] -1.75
## alpha_shelf[6] 2.02
## beta_sp[1] -3.57
## beta_sp[2] 2.98
## beta_sp[3] 2.98
## beta_sp[4] -11.17
## beta_sp[5] 6.62
## beta_sp_shade[1] 25.96
## beta_sp_shade[2] 14.97
## beta_sp_shade[3] 18.67
## beta_sp_shade[4] 25.55
## beta_sp_shade[5] 22.51
## sigma 13.04
## sigma_acs 7.53
## sigma_shelf 2.98
## beta_shade NA
## nobs 1008
## m.shade.species.int.acs.shelf
## alpha 41.19
## alpha_acs[1] -6.46
## alpha_acs[2] -11.28
## alpha_acs[3] -3.36
## alpha_acs[4] 3.41
## alpha_acs[5] 13.85
## alpha_acs[6] 2.73
## alpha_acs[7] 9.92
## alpha_acs[8] 6.94
## alpha_acs[9] 3.69
## alpha_acs[10] 2.99
## alpha_acs[11] 2.71
## alpha_acs[12] -10.11
## alpha_acs[13] -4.08
## alpha_acs[14] -0.78
## alpha_acs[15] -5.43
## alpha_acs[16] 1.86
## alpha_acs[17] 5.13
## alpha_acs[18] -12.23
## alpha_acs[19] 15.09
## alpha_acs[20] -5.10
## alpha_acs[21] 1.64
## alpha_acs[22] 4.55
## alpha_acs[23] -2.15
## alpha_acs[24] -4.07
## alpha_acs[25] -4.57
## alpha_acs[26] -0.08
## alpha_acs[27] 2.30
## alpha_acs[28] -4.6
## alpha_acs[29] 1.37
## alpha_acs[30] 6.84
## alpha_acs[31] 10.79
## alpha_acs[32] 1.29
## alpha_acs[33] -2.71
## alpha_acs[34] -8.68
## alpha_acs[35] -5.99
## alpha_acs[36] -6.42
## alpha_shelf[1] -3.09
## alpha_shelf[2] 3.62
## alpha_shelf[3] -0.15
## alpha_shelf[4] -2.37
## alpha_shelf[5] -1.05
## alpha_shelf[6] 2.70
## beta_sp[1] -3.72
## beta_sp[2] 2.74
## beta_sp[3] 2.81
## beta_sp[4] -11.21
## beta_sp[5] 6.58
## beta_sp_shade[1] 5.40
## beta_sp_shade[2] -5.33
## beta_sp_shade[3] -1.69
## beta_sp_shade[4] 5.05
## beta_sp_shade[5] 2.03
## sigma 13.05
## sigma_acs 7.54
## sigma_shelf 2.85
## beta_shade 22
## nobs 1008plot(coeftab(m.shade.species.sep.acs.shelf.no.alpha))
So long as an appropriate prior for the species intercepts is used this is OK. The other reasonable paramaterization would be to have the alpha befor the first species.
Consider different slopes for accessions
m.shade.species.sep.acs.int.shelf.no.alpha <- map2stan(
alist(
totleng ~ dnorm(mu,sigma),
mu <- alpha_acs[acs_id] +
alpha_shelf[shelf_id] +
beta_sp[species_id] +
beta_acs_shade[acs_id]*shade,
alpha_acs[acs_id] ~ dnorm(0,sigma_acs),
alpha_shelf[shelf_id] ~ dnorm(0,sigma_shelf),
beta_sp[species_id] ~ dnorm(42,10),
beta_acs_shade[acs_id] ~ dnorm(20,20),
c(sigma,sigma_acs,sigma_shelf) ~ dcauchy(0,1)
),
data = tomato.small,
iter = 4000, # Rhat a bit high with defaults
warmup = 1000,
chains=4,
cores=4
)##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 5e-06 seconds (Warm-up)
## 0.000264 seconds (Sampling)
## 0.000269 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 1200 / 12000 ]
[ 2400 / 12000 ]
[ 3600 / 12000 ]
[ 4800 / 12000 ]
[ 6000 / 12000 ]
[ 7200 / 12000 ]
[ 8400 / 12000 ]
[ 9600 / 12000 ]
[ 10800 / 12000 ]
[ 12000 / 12000 ]save.image("tomato1212.Rdata")plot(m.shade.species.sep.acs.int.shelf.no.alpha,ask=FALSE)## Waiting to draw page 2 of 6
## Waiting to draw page 3 of 6
## Waiting to draw page 4 of 6
## Waiting to draw page 5 of 6
## Waiting to draw page 6 of 6
par(mfrow=c(1,1))
precis(m.shade.species.sep.acs.int.shelf.no.alpha,depth = 2)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha_acs[1] -1.91 3.29 -7.36 3.14 12000 1
## alpha_acs[2] -4.09 3.23 -9.27 0.98 12000 1
## alpha_acs[3] -3.98 2.99 -8.62 0.86 12000 1
## alpha_acs[4] 2.81 3.03 -2.19 7.50 12000 1
## alpha_acs[5] 8.04 3.08 3.10 12.93 12000 1
## alpha_acs[6] -0.17 2.98 -4.98 4.42 12000 1
## alpha_acs[7] 5.83 3.19 0.39 10.57 12000 1
## alpha_acs[8] 0.10 2.83 -4.29 4.72 12000 1
## alpha_acs[9] 2.49 2.92 -1.93 7.39 12000 1
## alpha_acs[10] -0.01 3.41 -5.49 5.36 12000 1
## alpha_acs[11] 2.96 2.96 -1.85 7.58 12000 1
## alpha_acs[12] -4.95 3.08 -9.84 -0.04 12000 1
## alpha_acs[13] 0.24 3.25 -4.84 5.45 12000 1
## alpha_acs[14] 2.07 3.05 -2.62 7.16 12000 1
## alpha_acs[15] -1.00 3.99 -7.17 5.59 12000 1
## alpha_acs[16] 2.68 2.86 -1.84 7.27 12000 1
## alpha_acs[17] 4.75 3.01 -0.05 9.50 12000 1
## alpha_acs[18] -7.48 3.13 -12.23 -2.27 12000 1
## alpha_acs[19] 5.07 3.08 0.28 10.04 12000 1
## alpha_acs[20] -3.36 2.86 -7.96 1.15 12000 1
## alpha_acs[21] -0.10 2.81 -4.50 4.45 12000 1
## alpha_acs[22] 0.59 3.98 -5.74 6.80 12000 1
## alpha_acs[23] -1.73 3.40 -7.22 3.54 12000 1
## alpha_acs[24] -2.45 2.93 -7.04 2.29 12000 1
## alpha_acs[25] -1.85 2.87 -6.48 2.63 12000 1
## alpha_acs[26] -2.05 3.23 -7.31 2.92 12000 1
## alpha_acs[27] 1.06 2.84 -3.50 5.51 12000 1
## alpha_acs[28] -0.42 3.04 -5.24 4.41 12000 1
## alpha_acs[29] -0.57 3.09 -5.43 4.44 12000 1
## alpha_acs[30] 3.84 3.06 -1.16 8.57 12000 1
## alpha_acs[31] 5.43 3.10 0.58 10.37 12000 1
## alpha_acs[32] -0.17 3.15 -5.31 4.69 12000 1
## alpha_acs[33] 0.41 3.19 -4.55 5.59 12000 1
## alpha_acs[34] -6.79 3.13 -11.47 -1.49 12000 1
## alpha_acs[35] -4.06 3.15 -8.88 1.20 12000 1
## alpha_acs[36] -2.02 3.81 -8.13 3.92 12000 1
## alpha_shelf[1] -3.09 1.89 -6.13 -0.26 5976 1
## alpha_shelf[2] 4.21 1.92 1.17 7.17 5401 1
## alpha_shelf[3] 0.42 1.87 -2.44 3.45 5620 1
## alpha_shelf[4] -2.55 1.77 -5.34 0.28 6395 1
## alpha_shelf[5] -1.56 1.81 -4.37 1.32 6575 1
## alpha_shelf[6] 2.33 1.72 -0.43 4.98 6159 1
## beta_sp[1] 37.79 2.58 33.70 41.85 7218 1
## beta_sp[2] 44.56 2.55 40.53 48.58 7254 1
## beta_sp[3] 45.16 2.48 41.13 48.99 7411 1
## beta_sp[4] 29.76 2.75 25.41 34.14 8491 1
## beta_sp[5] 48.55 2.59 44.34 52.53 7240 1
## beta_acs_shade[1] 5.46 5.25 -3.05 13.81 12000 1
## beta_acs_shade[2] 9.11 4.81 1.54 16.90 12000 1
## beta_acs_shade[3] 17.70 4.37 10.65 24.55 12000 1
## beta_acs_shade[4] 15.83 4.43 8.78 22.93 12000 1
## beta_acs_shade[5] 25.37 4.39 18.43 32.48 8870 1
## beta_acs_shade[6] 21.90 4.59 14.47 29.14 12000 1
## beta_acs_shade[7] 29.12 4.57 21.92 36.47 12000 1
## beta_acs_shade[8] 33.91 4.27 27.29 40.95 12000 1
## beta_acs_shade[9] 19.82 4.34 12.87 26.63 12000 1
## beta_acs_shade[10] 34.55 5.71 25.29 43.41 12000 1
## beta_acs_shade[11] 21.08 4.30 13.97 27.54 12000 1
## beta_acs_shade[12] 7.11 4.93 -0.50 15.05 12000 1
## beta_acs_shade[13] 18.79 4.84 11.25 26.68 12000 1
## beta_acs_shade[14] 9.23 5.16 1.00 17.49 12000 1
## beta_acs_shade[15] 1.87 8.93 -12.13 16.42 12000 1
## beta_acs_shade[16] 19.99 4.14 13.56 26.80 12000 1
## beta_acs_shade[17] 26.30 4.45 19.32 33.55 12000 1
## beta_acs_shade[18] 17.92 4.76 9.99 25.11 12000 1
## beta_acs_shade[19] 41.74 4.54 34.26 48.73 12000 1
## beta_acs_shade[20] 14.91 4.30 8.19 21.88 12000 1
## beta_acs_shade[21] 22.02 4.40 15.35 29.35 12000 1
## beta_acs_shade[22] 35.38 6.96 24.00 46.06 12000 1
## beta_acs_shade[23] 27.05 5.08 18.78 34.89 12000 1
## beta_acs_shade[24] 12.63 4.31 5.52 19.31 12000 1
## beta_acs_shade[25] 10.03 4.15 3.02 16.26 12000 1
## beta_acs_shade[26] 31.52 5.00 23.88 39.86 12000 1
## beta_acs_shade[27] 29.19 4.28 22.35 35.93 12000 1
## beta_acs_shade[28] 17.81 4.59 10.53 25.20 12000 1
## beta_acs_shade[29] 31.20 4.81 23.50 38.93 12000 1
## beta_acs_shade[30] 31.58 4.49 24.27 38.49 12000 1
## beta_acs_shade[31] 27.52 4.73 19.87 34.93 12000 1
## beta_acs_shade[32] 31.18 4.95 23.01 38.84 12000 1
## beta_acs_shade[33] 21.13 4.61 13.78 28.37 12000 1
## beta_acs_shade[34] 20.46 4.59 12.94 27.48 12000 1
## beta_acs_shade[35] 19.27 5.06 10.90 27.17 12000 1
## beta_acs_shade[36] 19.25 6.42 8.04 28.47 12000 1
## sigma 12.72 0.29 12.25 13.19 12000 1
## sigma_acs 4.69 0.91 3.21 6.06 5566 1
## sigma_shelf 3.31 1.31 1.54 4.93 6730 1compare(m.shade.species.sep.acs.int.shelf.no.alpha,m.shade.species.sep.acs.shelf.no.alpha)## WAIC pWAIC dWAIC weight SE
## m.shade.species.sep.acs.int.shelf.no.alpha 8061.3 66.9 0.0 1 58.78
## m.shade.species.sep.acs.shelf.no.alpha 8083.5 42.4 22.3 0 56.47
## dSE
## m.shade.species.sep.acs.int.shelf.no.alpha NA
## m.shade.species.sep.acs.shelf.no.alpha 20.4coeftab(m.shade.species.sep.acs.int.shelf.no.alpha,m.shade.species.sep.acs.shelf.no.alpha)## m.shade.species.sep.acs.int.shelf.no.alpha
## alpha_acs[1] -1.91
## alpha_acs[2] -4.09
## alpha_acs[3] -3.98
## alpha_acs[4] 2.81
## alpha_acs[5] 8.04
## alpha_acs[6] -0.17
## alpha_acs[7] 5.83
## alpha_acs[8] 0.10
## alpha_acs[9] 2.49
## alpha_acs[10] -0.01
## alpha_acs[11] 2.96
## alpha_acs[12] -4.95
## alpha_acs[13] 0.24
## alpha_acs[14] 2.07
## alpha_acs[15] -1.00
## alpha_acs[16] 2.68
## alpha_acs[17] 4.75
## alpha_acs[18] -7.48
## alpha_acs[19] 5.07
## alpha_acs[20] -3.36
## alpha_acs[21] -0.10
## alpha_acs[22] 0.59
## alpha_acs[23] -1.73
## alpha_acs[24] -2.45
## alpha_acs[25] -1.85
## alpha_acs[26] -2.05
## alpha_acs[27] 1.06
## alpha_acs[28] -0.42
## alpha_acs[29] -0.57
## alpha_acs[30] 3.84
## alpha_acs[31] 5.43
## alpha_acs[32] -0.17
## alpha_acs[33] 0.41
## alpha_acs[34] -6.79
## alpha_acs[35] -4.06
## alpha_acs[36] -2.02
## alpha_shelf[1] -3.09
## alpha_shelf[2] 4.21
## alpha_shelf[3] 0.42
## alpha_shelf[4] -2.55
## alpha_shelf[5] -1.56
## alpha_shelf[6] 2.33
## beta_sp[1] 37.79
## beta_sp[2] 44.56
## beta_sp[3] 45.16
## beta_sp[4] 29.76
## beta_sp[5] 48.55
## beta_acs_shade[1] 5.46
## beta_acs_shade[2] 9.11
## beta_acs_shade[3] 17.7
## beta_acs_shade[4] 15.83
## beta_acs_shade[5] 25.37
## beta_acs_shade[6] 21.9
## beta_acs_shade[7] 29.12
## beta_acs_shade[8] 33.91
## beta_acs_shade[9] 19.82
## beta_acs_shade[10] 34.55
## beta_acs_shade[11] 21.08
## beta_acs_shade[12] 7.11
## beta_acs_shade[13] 18.79
## beta_acs_shade[14] 9.23
## beta_acs_shade[15] 1.87
## beta_acs_shade[16] 19.99
## beta_acs_shade[17] 26.3
## beta_acs_shade[18] 17.92
## beta_acs_shade[19] 41.74
## beta_acs_shade[20] 14.91
## beta_acs_shade[21] 22.02
## beta_acs_shade[22] 35.38
## beta_acs_shade[23] 27.05
## beta_acs_shade[24] 12.63
## beta_acs_shade[25] 10.03
## beta_acs_shade[26] 31.52
## beta_acs_shade[27] 29.19
## beta_acs_shade[28] 17.81
## beta_acs_shade[29] 31.2
## beta_acs_shade[30] 31.58
## beta_acs_shade[31] 27.52
## beta_acs_shade[32] 31.18
## beta_acs_shade[33] 21.13
## beta_acs_shade[34] 20.46
## beta_acs_shade[35] 19.27
## beta_acs_shade[36] 19.25
## sigma 12.72
## sigma_acs 4.69
## sigma_shelf 3.31
## beta_sp_shade[1] NA
## beta_sp_shade[2] NA
## beta_sp_shade[3] NA
## beta_sp_shade[4] NA
## beta_sp_shade[5] NA
## nobs 1008
## m.shade.species.sep.acs.shelf.no.alpha
## alpha_acs[1] -6.56
## alpha_acs[2] -11.70
## alpha_acs[3] -3.38
## alpha_acs[4] 3.34
## alpha_acs[5] 13.82
## alpha_acs[6] 2.64
## alpha_acs[7] 9.60
## alpha_acs[8] 6.77
## alpha_acs[9] 3.54
## alpha_acs[10] 3.60
## alpha_acs[11] 2.32
## alpha_acs[12] -10.32
## alpha_acs[13] -3.39
## alpha_acs[14] -0.94
## alpha_acs[15] -5.62
## alpha_acs[16] 1.46
## alpha_acs[17] 5.19
## alpha_acs[18] -12.21
## alpha_acs[19] 14.72
## alpha_acs[20] -5.28
## alpha_acs[21] 1.46
## alpha_acs[22] 5.18
## alpha_acs[23] -1.51
## alpha_acs[24] -4.11
## alpha_acs[25] -4.64
## alpha_acs[26] -0.02
## alpha_acs[27] 2.34
## alpha_acs[28] -4.56
## alpha_acs[29] 1.42
## alpha_acs[30] 6.90
## alpha_acs[31] 10.64
## alpha_acs[32] 1.98
## alpha_acs[33] -1.99
## alpha_acs[34] -9.05
## alpha_acs[35] -6.39
## alpha_acs[36] -5.79
## alpha_shelf[1] -2.40
## alpha_shelf[2] 4.37
## alpha_shelf[3] 0.57
## alpha_shelf[4] -3.03
## alpha_shelf[5] -1.68
## alpha_shelf[6] 2.11
## beta_sp[1] 38.01
## beta_sp[2] 44.79
## beta_sp[3] 44.91
## beta_sp[4] 29.75
## beta_sp[5] 48.92
## beta_acs_shade[1] NA
## beta_acs_shade[2] NA
## beta_acs_shade[3] NA
## beta_acs_shade[4] NA
## beta_acs_shade[5] NA
## beta_acs_shade[6] NA
## beta_acs_shade[7] NA
## beta_acs_shade[8] NA
## beta_acs_shade[9] NA
## beta_acs_shade[10] NA
## beta_acs_shade[11] NA
## beta_acs_shade[12] NA
## beta_acs_shade[13] NA
## beta_acs_shade[14] NA
## beta_acs_shade[15] NA
## beta_acs_shade[16] NA
## beta_acs_shade[17] NA
## beta_acs_shade[18] NA
## beta_acs_shade[19] NA
## beta_acs_shade[20] NA
## beta_acs_shade[21] NA
## beta_acs_shade[22] NA
## beta_acs_shade[23] NA
## beta_acs_shade[24] NA
## beta_acs_shade[25] NA
## beta_acs_shade[26] NA
## beta_acs_shade[27] NA
## beta_acs_shade[28] NA
## beta_acs_shade[29] NA
## beta_acs_shade[30] NA
## beta_acs_shade[31] NA
## beta_acs_shade[32] NA
## beta_acs_shade[33] NA
## beta_acs_shade[34] NA
## beta_acs_shade[35] NA
## beta_acs_shade[36] NA
## sigma 13.04
## sigma_acs 7.59
## sigma_shelf 2.97
## beta_sp_shade[1] 26.07
## beta_sp_shade[2] 15.01
## beta_sp_shade[3] 18.76
## beta_sp_shade[4] 25.85
## beta_sp_shade[5] 22.48
## nobs 1008plot(coeftab(m.shade.species.sep.acs.int.shelf.no.alpha,m.shade.species.sep.acs.shelf.no.alpha),pars=c(1,44:56))
drop the species shade interaction
m.shade.species.acs.int.shelf <- map2stan(
alist(
totleng ~ dnorm(mu,sigma),
mu <- alpha +
alpha_acs[acs_id] +
alpha_shelf[shelf_id] +
beta_sp[species_id] +
beta_shade*shade +
beta_acs_shade[acs_id]*shade,
alpha ~ dnorm(53,20),
alpha_acs[acs_id] ~ dnorm(0,sigma_acs),
alpha_shelf[shelf_id] ~ dnorm(0,sigma_shelf),
beta_sp[species_id] ~ dnorm(0,10),
beta_shade ~ dnorm(0,20),
beta_acs_shade[acs_id] ~ dnorm(0,10),
c(sigma,sigma_acs,sigma_shelf) ~ dcauchy(0,1)
),
data = tomato.small,
iter = 4000, # Rhat a bit high with defaults
warmup = 1000,
chains=4,
cores=4
)##
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 5e-06 seconds (Warm-up)
## 0.00032 seconds (Sampling)
## 0.000325 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 1200 / 12000 ]
[ 2400 / 12000 ]
[ 3600 / 12000 ]
[ 4800 / 12000 ]
[ 6000 / 12000 ]
[ 7200 / 12000 ]
[ 8400 / 12000 ]
[ 9600 / 12000 ]
[ 10800 / 12000 ]
[ 12000 / 12000 ]save.image("tomato1212.Rdata")plot(m.shade.species.acs.int.shelf,ask=FALSE)## Waiting to draw page 2 of 6
## Waiting to draw page 3 of 6
## Waiting to draw page 4 of 6
## Waiting to draw page 5 of 6
## Waiting to draw page 6 of 6
par(mfrow=c(1,1))
precis(m.shade.species.acs.int.shelf,depth = 2)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha 41.63 5.03 33.76 49.41 2006 1.00
## alpha_acs[1] -2.66 3.37 -8.13 2.46 12000 1.00
## alpha_acs[2] -5.03 3.31 -10.20 0.23 12000 1.00
## alpha_acs[3] -4.04 3.11 -8.93 0.94 12000 1.00
## alpha_acs[4] 2.96 3.11 -1.94 7.91 12000 1.00
## alpha_acs[5] 8.81 3.21 3.55 13.73 12000 1.00
## alpha_acs[6] 0.16 3.07 -4.96 4.85 12000 1.00
## alpha_acs[7] 6.44 3.23 1.23 11.51 12000 1.00
## alpha_acs[8] 0.93 2.90 -3.57 5.61 12000 1.00
## alpha_acs[9] 2.86 2.99 -1.92 7.70 12000 1.00
## alpha_acs[10] 0.36 3.49 -5.02 6.03 12000 1.00
## alpha_acs[11] 3.07 3.05 -1.85 7.87 12000 1.00
## alpha_acs[12] -5.59 3.13 -10.46 -0.48 12000 1.00
## alpha_acs[13] -0.25 3.37 -5.61 5.14 12000 1.00
## alpha_acs[14] 1.85 3.12 -3.00 6.87 12000 1.00
## alpha_acs[15] -2.25 4.07 -8.42 4.46 12000 1.00
## alpha_acs[16] 2.71 2.99 -2.13 7.31 12000 1.00
## alpha_acs[17] 4.89 3.08 -0.23 9.54 12000 1.00
## alpha_acs[18] -8.25 3.24 -13.46 -3.12 12000 1.00
## alpha_acs[19] 6.31 3.20 1.35 11.48 12000 1.00
## alpha_acs[20] -3.39 2.94 -7.89 1.44 12000 1.00
## alpha_acs[21] 0.27 2.86 -4.13 4.92 12000 1.00
## alpha_acs[22] 1.56 4.09 -4.95 7.93 12000 1.00
## alpha_acs[23] -1.86 3.52 -7.53 3.60 12000 1.00
## alpha_acs[24] -2.66 3.04 -7.78 1.91 12000 1.00
## alpha_acs[25] -2.06 3.01 -6.89 2.70 12000 1.00
## alpha_acs[26] -1.76 3.37 -7.01 3.68 12000 1.00
## alpha_acs[27] 1.15 3.03 -3.61 6.02 12000 1.00
## alpha_acs[28] -0.85 3.13 -5.94 3.99 12000 1.00
## alpha_acs[29] -0.30 3.15 -5.36 4.70 12000 1.00
## alpha_acs[30] 4.29 3.17 -0.98 9.13 12000 1.00
## alpha_acs[31] 6.32 3.29 1.04 11.56 12000 1.00
## alpha_acs[32] -0.06 3.23 -5.15 5.10 12000 1.00
## alpha_acs[33] 0.01 3.31 -5.26 5.28 12000 1.00
## alpha_acs[34] -7.19 3.17 -12.24 -2.17 12000 1.00
## alpha_acs[35] -4.38 3.24 -9.72 0.63 12000 1.00
## alpha_acs[36] -2.62 3.98 -8.63 3.97 12000 1.00
## alpha_shelf[1] -3.32 2.65 -7.02 0.18 2576 1.00
## alpha_shelf[2] 4.03 2.67 0.36 7.61 2518 1.00
## alpha_shelf[3] 0.14 2.62 -3.33 3.85 2538 1.00
## alpha_shelf[4] -2.50 2.67 -6.09 1.09 705 1.01
## alpha_shelf[5] -1.45 2.70 -5.05 2.15 733 1.00
## alpha_shelf[6] 2.44 2.62 -0.76 6.27 763 1.01
## beta_sp[1] -3.67 4.79 -11.25 4.02 5127 1.00
## beta_sp[2] 2.38 4.77 -4.92 10.10 4957 1.00
## beta_sp[3] 2.96 4.72 -4.57 10.43 5007 1.00
## beta_sp[4] -11.43 4.82 -19.00 -3.61 5223 1.00
## beta_sp[5] 6.82 4.77 -0.89 14.26 4899 1.00
## beta_shade 22.00 3.91 16.43 27.60 1484 1.00
## beta_acs_shade[1] -13.41 4.79 -20.71 -5.58 12000 1.00
## beta_acs_shade[2] -10.63 4.43 -17.62 -3.38 12000 1.00
## beta_acs_shade[3] -3.22 4.16 -9.60 3.70 12000 1.00
## beta_acs_shade[4] -5.11 4.08 -11.72 1.31 12000 1.00
## beta_acs_shade[5] 3.25 4.18 -3.49 9.80 12000 1.00
## beta_acs_shade[6] 0.41 4.21 -6.09 7.27 12000 1.00
## beta_acs_shade[7] 6.49 4.22 -0.15 13.29 12000 1.00
## beta_acs_shade[8] 11.12 4.00 4.74 17.49 12000 1.00
## beta_acs_shade[9] -1.59 4.13 -7.99 5.08 12000 1.00
## beta_acs_shade[10] 10.61 5.12 2.69 19.03 12000 1.00
## beta_acs_shade[11] -0.63 4.03 -7.08 5.75 12000 1.00
## beta_acs_shade[12] -12.23 4.49 -19.65 -5.26 12000 1.00
## beta_acs_shade[13] -2.65 4.47 -9.67 4.58 12000 1.00
## beta_acs_shade[14] -10.53 4.72 -17.92 -2.87 12000 1.00
## beta_acs_shade[15] -12.12 7.13 -23.57 -0.85 12000 1.00
## beta_acs_shade[16] -1.54 3.87 -7.84 4.49 12000 1.00
## beta_acs_shade[17] 3.93 4.12 -2.53 10.54 12000 1.00
## beta_acs_shade[18] -2.97 4.45 -10.05 4.23 12000 1.00
## beta_acs_shade[19] 17.59 4.22 11.00 24.49 12000 1.00
## beta_acs_shade[20] -5.83 4.02 -12.27 0.58 12000 1.00
## beta_acs_shade[21] 0.49 4.05 -6.20 6.69 12000 1.00
## beta_acs_shade[22] 10.19 6.04 0.72 19.80 12000 1.00
## beta_acs_shade[23] 4.62 4.74 -2.94 12.18 12000 1.00
## beta_acs_shade[24] -7.81 4.03 -14.20 -1.41 12000 1.00
## beta_acs_shade[25] -10.32 3.91 -16.49 -4.10 12000 1.00
## beta_acs_shade[26] 8.58 4.62 1.31 16.05 12000 1.00
## beta_acs_shade[27] 6.74 4.02 0.44 13.24 12000 1.00
## beta_acs_shade[28] -3.42 4.30 -10.48 3.22 12000 1.00
## beta_acs_shade[29] 8.32 4.37 1.33 15.24 12000 1.00
## beta_acs_shade[30] 8.68 4.13 1.82 15.00 12000 1.00
## beta_acs_shade[31] 5.15 4.45 -1.93 12.30 12000 1.00
## beta_acs_shade[32] 8.08 4.56 1.02 15.55 12000 1.00
## beta_acs_shade[33] -0.59 4.31 -7.20 6.55 12000 1.00
## beta_acs_shade[34] -0.76 4.26 -7.54 5.99 12000 1.00
## beta_acs_shade[35] -1.88 4.67 -9.39 5.46 12000 1.00
## beta_acs_shade[36] -1.88 5.71 -10.77 7.41 12000 1.00
## sigma 12.72 0.29 12.24 13.17 12000 1.00
## sigma_acs 5.06 0.95 3.62 6.55 6622 1.00
## sigma_shelf 3.64 1.95 1.49 5.49 620 1.01compare(m.shade.species.int.acs.shelf,m.shade.species.acs.int.shelf)## WAIC pWAIC dWAIC weight SE dSE
## m.shade.species.acs.int.shelf 8056.3 63.7 0 1 58.40 NA
## m.shade.species.int.acs.shelf 8083.3 42.2 27 0 56.41 17.29coeftab(m.shade.species.acs.int.shelf)## m.shade.species.acs.int.shelf
## alpha 41.63
## alpha_acs[1] -2.66
## alpha_acs[2] -5.03
## alpha_acs[3] -4.04
## alpha_acs[4] 2.96
## alpha_acs[5] 8.81
## alpha_acs[6] 0.16
## alpha_acs[7] 6.44
## alpha_acs[8] 0.93
## alpha_acs[9] 2.86
## alpha_acs[10] 0.36
## alpha_acs[11] 3.07
## alpha_acs[12] -5.59
## alpha_acs[13] -0.25
## alpha_acs[14] 1.85
## alpha_acs[15] -2.25
## alpha_acs[16] 2.71
## alpha_acs[17] 4.89
## alpha_acs[18] -8.25
## alpha_acs[19] 6.31
## alpha_acs[20] -3.39
## alpha_acs[21] 0.27
## alpha_acs[22] 1.56
## alpha_acs[23] -1.86
## alpha_acs[24] -2.66
## alpha_acs[25] -2.06
## alpha_acs[26] -1.76
## alpha_acs[27] 1.15
## alpha_acs[28] -0.85
## alpha_acs[29] -0.3
## alpha_acs[30] 4.29
## alpha_acs[31] 6.32
## alpha_acs[32] -0.06
## alpha_acs[33] 0.01
## alpha_acs[34] -7.19
## alpha_acs[35] -4.38
## alpha_acs[36] -2.62
## alpha_shelf[1] -3.32
## alpha_shelf[2] 4.03
## alpha_shelf[3] 0.14
## alpha_shelf[4] -2.5
## alpha_shelf[5] -1.45
## alpha_shelf[6] 2.44
## beta_sp[1] -3.67
## beta_sp[2] 2.38
## beta_sp[3] 2.96
## beta_sp[4] -11.43
## beta_sp[5] 6.82
## beta_shade 22
## beta_acs_shade[1] -13.41
## beta_acs_shade[2] -10.63
## beta_acs_shade[3] -3.22
## beta_acs_shade[4] -5.11
## beta_acs_shade[5] 3.25
## beta_acs_shade[6] 0.41
## beta_acs_shade[7] 6.49
## beta_acs_shade[8] 11.12
## beta_acs_shade[9] -1.59
## beta_acs_shade[10] 10.61
## beta_acs_shade[11] -0.63
## beta_acs_shade[12] -12.23
## beta_acs_shade[13] -2.65
## beta_acs_shade[14] -10.53
## beta_acs_shade[15] -12.12
## beta_acs_shade[16] -1.54
## beta_acs_shade[17] 3.93
## beta_acs_shade[18] -2.97
## beta_acs_shade[19] 17.59
## beta_acs_shade[20] -5.83
## beta_acs_shade[21] 0.49
## beta_acs_shade[22] 10.19
## beta_acs_shade[23] 4.62
## beta_acs_shade[24] -7.81
## beta_acs_shade[25] -10.32
## beta_acs_shade[26] 8.58
## beta_acs_shade[27] 6.74
## beta_acs_shade[28] -3.42
## beta_acs_shade[29] 8.32
## beta_acs_shade[30] 8.68
## beta_acs_shade[31] 5.15
## beta_acs_shade[32] 8.08
## beta_acs_shade[33] -0.59
## beta_acs_shade[34] -0.76
## beta_acs_shade[35] -1.88
## beta_acs_shade[36] -1.88
## sigma 12.72
## sigma_acs 5.06
## sigma_shelf 3.64
## nobs 1008plot(coeftab(m.shade.species.acs.int.shelf),pars=c(1:42))
plot(coeftab(m.shade.species.acs.int.shelf),pars=c(43:86))
We don’t need both a species response and an accession response (interaction) term. If we are only keeping one of the terms then the accession interaction model is preferred.
From this I conclude that variation is better considered at the accession level than at the species level.
Smoking deaths among doctors
In 1961 Doll and Hill sent out a questionnaire to all men on the British Medical Register inquiring about their smoking habits. Almost 70% of such men replied. Death certificates were obtained for medical practitioners and causes of death were assigned on the basis of these certificates. The breslow data set contains the person-years of observations and deaths from coronary artery disease accumulated during the first ten years of the study.
Analyse this data set to determine the posterior probability that smoking increases death by coronary artery disease, that age increases death by coronary artery disease, and that there is an interaction between age and smoking.
You can load the data set and learn about the columns using the commands below
data("breslow",package = "boot")
help("breslow",package ="boot")
breslow## age smoke n y ns
## 1 40 0 18790 2 0
## 2 50 0 10673 12 0
## 3 60 0 5710 28 0
## 4 70 0 2585 28 0
## 5 80 0 1462 31 0
## 6 40 1 52407 32 52407
## 7 50 1 43248 104 43248
## 8 60 1 28612 206 28612
## 9 70 1 12663 186 12663
## 10 80 1 5317 102 5317You can think of “person-years” as the number of observations
breslow$n <- as.integer(breslow$n)
m.smoke <- map2stan(alist(
y~dbinom(n,p),
logit(p) <- alpha + beta_smoke * smoke,
alpha ~ dnorm(0,5),
beta_smoke ~ dnorm(0,5)),
data=breslow,
chains = 4,
cores = 2)## Warning in FUN(X[[i]], ...): data with name age is not numeric and not used
## Warning in FUN(X[[i]], ...): data with name age is not numeric and not used##
## SAMPLING FOR MODEL 'y ~ dbinom(n, p)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 3e-06 seconds (Warm-up)
## 3.2e-05 seconds (Sampling)
## 3.5e-05 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.save.image("tomato1212.Rdata")plot(m.smoke)
pairs(m.smoke)
precis(m.smoke)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha -5.96 0.10 -6.12 -5.81 843 1.01
## beta_smoke 0.55 0.11 0.39 0.73 844 1.01logistic(-5.97)## [1] 0.002547734logistic(-5.97 + 0.55)## [1] 0.004407633breslow$age_id <- coerce_index(breslow$age)
m.age <- map2stan(alist(
y~dbinom(n,p),
logit(p) <- alpha_age[age_id],
alpha_age[age_id] ~ dnorm(0,5)),
data=breslow,
chains = 4,
cores = 2)## Warning in FUN(X[[i]], ...): data with name age is not numeric and not used
## Warning in FUN(X[[i]], ...): data with name age is not numeric and not used##
## SAMPLING FOR MODEL 'y ~ dbinom(n, p)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 5e-06 seconds (Warm-up)
## 3.7e-05 seconds (Sampling)
## 4.2e-05 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.save.image("tomato1212.Rdata")plot(m.age)
precis(m.age,depth=2)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha_age[1] -7.65 0.17 -7.91 -7.38 3918 1
## alpha_age[2] -6.14 0.09 -6.28 -6.00 4000 1
## alpha_age[3] -4.98 0.06 -5.09 -4.88 4000 1
## alpha_age[4] -4.25 0.07 -4.36 -4.15 3259 1
## alpha_age[5] -3.91 0.08 -4.05 -3.78 4000 1logistic(precis(m.age,depth=2)@output$Mean)## [1] 0.0004743055 0.0021526771 0.0068113465 0.0140058684 0.0195927144pairs(m.age)
compare(m.age,m.smoke)## WAIC pWAIC dWAIC weight SE dSE
## m.age 8614.6 4.6 0.0 1 268.90 NA
## m.smoke 9495.9 2.0 881.3 0 296.64 62.33m.smoke.age <- map2stan(
alist(y~dbinom(n,p),
logit(p) <- alpha_age[age_id] + beta_smoke*smoke,
alpha_age[age_id] ~ dnorm(0,5),
beta_smoke ~ dnorm(0,5)),
data=breslow,
chains = 4,
cores = 2)## Warning in FUN(X[[i]], ...): data with name age is not numeric and not used
## Warning in FUN(X[[i]], ...): data with name age is not numeric and not used##
## SAMPLING FOR MODEL 'y ~ dbinom(n, p)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 4e-06 seconds (Warm-up)
## 3.7e-05 seconds (Sampling)
## 4.1e-05 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.save.image("tomato1212.Rdata")plot(m.smoke.age)
pairs(m.smoke.age)
precis(m.smoke.age,depth=2)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha_age[1] -7.92 0.19 -8.20 -7.61 1590 1.00
## alpha_age[2] -6.43 0.13 -6.63 -6.23 1133 1.00
## alpha_age[3] -5.28 0.11 -5.45 -5.10 1018 1.01
## alpha_age[4] -4.55 0.11 -4.73 -4.38 1095 1.00
## alpha_age[5] -4.20 0.12 -4.39 -4.00 1098 1.01
## beta_smoke 0.35 0.10 0.17 0.50 839 1.01compare(m.age,m.smoke,m.smoke.age)## WAIC pWAIC dWAIC weight SE dSE
## m.smoke.age 8604.9 5.7 0.0 0.99 268.48 NA
## m.age 8614.6 4.6 9.7 0.01 268.90 6.56
## m.smoke 9495.9 2.0 891.0 0.00 296.64 61.11m.smoke.age.int <- map2stan(
alist(y~dbinom(n,p),
logit(p) <- alpha_age[age_id] +
beta_smoke_age[age_id]*smoke,
alpha_age[age_id] ~ dnorm(0,5),
beta_smoke_age[age_id] ~ dnorm(0,5)),
data=breslow,
chains = 4,
cores = 2)## Warning in FUN(X[[i]], ...): data with name age is not numeric and not used
## Warning in FUN(X[[i]], ...): data with name age is not numeric and not used##
## SAMPLING FOR MODEL 'y ~ dbinom(n, p)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
<<<<<<< HEAD
## Elapsed Time: 8e-06 seconds (Warm-up)
## 3.9e-05 seconds (Sampling)
## 4.7e-05 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.save.image("tomato1212.Rdata")plot(m.smoke.age.int)
precis(m.smoke.age.int,depth=2)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
<<<<<<< HEAD
## alpha_age[1] -9.16 0.65 -10.12 -8.16 2394 1
## alpha_age[2] -6.80 0.29 -7.23 -6.32 3034 1
## alpha_age[3] -5.32 0.18 -5.58 -5.01 3019 1
## alpha_age[4] -4.53 0.18 -4.80 -4.24 2993 1
## alpha_age[5] -3.85 0.18 -4.13 -3.57 3198 1
## beta_smoke_age[1] 1.74 0.67 0.70 2.77 2493 1
## beta_smoke_age[2] 0.77 0.30 0.29 1.26 3083 1
## beta_smoke_age[3] 0.39 0.19 0.08 0.69 3092 1
## beta_smoke_age[4] 0.32 0.19 0.02 0.62 2762 1
## beta_smoke_age[5] -0.09 0.20 -0.41 0.23 3144 1compare(m.smoke,m.age,m.smoke.age,m.smoke.age.int)## WAIC pWAIC dWAIC weight SE dSE
## m.smoke.age.int 8600.5 9.5 0.0 0.89 268.30 NA
## m.smoke.age 8604.7 5.6 4.2 0.11 268.47 6.83
## m.age 8614.5 4.5 14.0 0.00 268.92 9.32
## m.smoke 9495.8 1.9 895.3 0.00 296.55 61.26
compare(m.smoke,m.age,m.smoke.age,m.smoke.age.int)## WAIC pWAIC dWAIC weight SE dSE
## m.smoke.age.int 8600.8 9.7 0.0 0.88 268.39 NA
## m.smoke.age 8604.9 5.7 4.1 0.11 268.48 6.94
## m.age 8614.6 4.6 13.8 0.00 268.90 9.36
## m.smoke 9495.9 2.0 895.1 0.00 296.64 61.33
pred <- ensemble(m.smoke.age,m.smoke.age.int)
mu <- apply(pred$link,2,mean)
PI <- apply(pred$link,2,PI)
plot.frame <- data.frame(age=breslow$age,smoke=as.factor(breslow$smoke),mean=mu,low=PI[1,],high=PI[2,])
plot.frame## age smoke mean low high
<<<<<<< HEAD
## 1 40 0 0.0001500566 3.539396e-05 0.0003566858
## 2 50 0 0.0011932146 6.774274e-04 0.0017858632
## 3 60 0 0.0049543903 3.790866e-03 0.0064157975
## 4 70 0 0.0108018511 8.230865e-03 0.0138126302
## 5 80 0 0.0205305623 1.406183e-02 0.0271009650
## 6 40 1 0.0005962079 4.365823e-04 0.0007722487
## 7 50 1 0.0023991355 2.032179e-03 0.0027737890
## 8 60 1 0.0071879859 6.397410e-03 0.0080402895
## 9 70 1 0.0147819631 1.307582e-02 0.0165350039
## 10 80 1 0.0193259024 1.655181e-02 0.0222723804pl <- ggplot(plot.frame,aes(x=age,y=mean,ymax=high,ymin=low,fill=smoke))
pl <- pl + geom_bar(position = "dodge",stat = "identity")
pl + geom_errorbar(position = position_dodge(width=0.9),width=0.5)
pl <- ggplot(plot.frame,aes(x=age,y=mean,ymax=high,ymin=low,fill=smoke))
pl <- pl + geom_bar(position = "dodge",stat = "identity")
pl + geom_errorbar(position = position_dodge(width=0.9),width=0.5)
Predict from interaction model only:
pred <- link(m.smoke.age.int)## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]mu <- apply(pred,2,mean)
PI <- apply(pred,2,PI)
plot.frame <- data.frame(age=breslow$age,smoke=as.factor(breslow$smoke),mean=mu,low=PI[1,],high=PI[2,])
plot.frame## age smoke mean low high
## 1 40 0 0.0001256134 0.0000304405 0.0002708942
## 2 50 0 0.0011638369 0.0006797473 0.0017323160
## 3 60 0 0.0049193605 0.0036554407 0.0064815203
## 4 70 0 0.0109677605 0.0079951939 0.0144062696
## 5 80 0 0.0214162370 0.0157445032 0.0272010499
## 6 40 1 0.0006083455 0.0004577301 0.0007929347
## 7 50 1 0.0023939922 0.0020190133 0.0027827213
## 8 60 1 0.0071691387 0.0064104530 0.0079867875
## 9 70 1 0.0146927590 0.0129661484 0.0164295508
## 10 80 1 0.0191886953 0.0163124232 0.0223023895pl <- ggplot(plot.frame,aes(x=age,y=mean,ymax=high,ymin=low,fill=smoke))
pl <- pl + geom_bar(position = "dodge",stat = "identity")
pl + geom_errorbar(position = position_dodge(width=0.9),width=0.5)
breslow$age_num <- as.numeric(as.character(breslow$age))
m.smoke.age.int.numeric <- map2stan(
alist(y~dbinom(n,p),
logit(p) <- alpha +
beta_age * age_num +
beta_smoke*smoke +
beta_smoke_age*smoke*age_num,
alpha ~ dnorm(0,10),
c(beta_smoke,beta_age) ~ dnorm(0,5),
beta_smoke_age ~ dnorm(0,1)),
data=breslow,
chains = 4,
cores = 2)## Warning in FUN(X[[i]], ...): data with name age is not numeric and not used
## Warning in FUN(X[[i]], ...): data with name age is not numeric and not used##
## SAMPLING FOR MODEL 'y ~ dbinom(n, p)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 3e-06 seconds (Warm-up)
## 3.3e-05 seconds (Sampling)
## 3.6e-05 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.save.image("tomato1212.Rdata")plot(m.smoke.age.int.numeric)
precis(m.smoke.age.int.numeric)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha -12.07 0.52 -12.95 -11.29 270 1.01
## beta_smoke 2.07 0.56 1.19 2.96 266 1.01
## beta_age 0.11 0.01 0.09 0.12 271 1.01
## beta_smoke_age -0.03 0.01 -0.04 -0.01 268 1.01compare(m.smoke.age.int, m.smoke.age.int.numeric)## WAIC pWAIC dWAIC weight SE dSE
## m.smoke.age.int 8600.5 9.5 0.0 1 268.30 NA
## m.smoke.age.int.numeric 8647.8 3.5 47.3 0 269.84 15
## WAIC pWAIC dWAIC weight SE dSE
## m.smoke.age.int 8600.8 9.7 0.0 1 268.39 NA
## m.smoke.age.int.numeric 8647.4 3.3 46.5 0 269.81 14.99
Cane Sugar
This data comes from an experiment to measure disease resistance in different varieties of sugar cane.
You can get the data and learn about it with:
data("cane",package="boot")
help("cane",package="boot")
head(cane)## n r x var block
## 1 87 76 19 1 A
## 2 119 8 14 2 A
## 3 94 74 9 3 A
## 4 95 11 12 4 A
## 5 134 0 12 5 A
## 6 92 0 3 6 Asummary(cane)## n r x var block
## Min. : 29.00 Min. : 0.00 Min. : 1.00 1 : 4 A:45
## 1st Qu.: 85.75 1st Qu.: 3.00 1st Qu.: 7.00 10 : 4 B:45
## Median :112.50 Median : 11.50 Median :11.50 11 : 4 C:45
## Mean :118.14 Mean : 20.26 Mean :11.94 12 : 4 D:45
## 3rd Qu.:144.50 3rd Qu.: 25.25 3rd Qu.:15.00 13 : 4
## Max. :243.00 Max. :131.00 Max. :36.00 14 : 4
## (Other):156Is there evidence of differences in disease resistance in the different varieties? Does including an adaptive prior for Block improve the model fit?
mean(cane$r/cane$n)## [1] 0.1741397logit(mean(cane$r/cane$n))## [1] -1.556568cane$var_id <- coerce_index(cane$var)
cane$block_id <- coerce_index(cane$block)
cane$n <- as.integer(cane$n)
m.cane.alpha <- map2stan(
alist(
r ~ dbinom(n,p),
logit(p) <- alpha,
alpha <- dnorm(-1.55,1)
),
data=cane,
chains = 4,
cores = 2
)## Warning in FUN(X[[i]], ...): data with name var is not numeric and not used## Warning in FUN(X[[i]], ...): data with name block is not numeric and not
## used## Warning in FUN(X[[i]], ...): data with name var is not numeric and not used## Warning in FUN(X[[i]], ...): data with name block is not numeric and not
## used##
## SAMPLING FOR MODEL 'r ~ dbinom(n, p)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 3e-06 seconds (Warm-up)
## 0.00014 seconds (Sampling)
## 0.000143 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.plot(m.cane.alpha)
precis(m.cane.alpha)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha -1.58 0.02 -1.6 -1.55 3284 1logistic(-1.58)## [1] 0.1707955
m.var <- map2stan(
alist(
r ~ dbinom(n,p),
logit(p) <- alpha + alpha_var[var_id],
alpha <- dnorm(-1.55,1),
alpha_var[var_id] ~ dnorm(0,5)
),
data=cane,
chains = 4,
cores = 2,
iter=10000,
warmup = 1000
)## Warning in FUN(X[[i]], ...): data with name var is not numeric and not used## Warning in FUN(X[[i]], ...): data with name block is not numeric and not
## used## Warning in FUN(X[[i]], ...): data with name var is not numeric and not used## Warning in FUN(X[[i]], ...): data with name block is not numeric and not
## used##
## SAMPLING FOR MODEL 'r ~ dbinom(n, p)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 3e-06 seconds (Warm-up)
## 0.000159 seconds (Sampling)
## 0.000162 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 3600 / 36000 ]
[ 7200 / 36000 ]
[ 10800 / 36000 ]
[ 14400 / 36000 ]
[ 18000 / 36000 ]
[ 21600 / 36000 ]
[ 25200 / 36000 ]
[ 28800 / 36000 ]
[ 32400 / 36000 ]
[ 36000 / 36000 ]## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.save.image("tomato1212.Rdata")plot(m.var,ask=FALSE)## Waiting to draw page 2 of 4
## Waiting to draw page 3 of 4
## Waiting to draw page 4 of 4
precis(m.var,depth=2)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha -1.93 0.60 -2.90 -0.96 953 1.01
## alpha_var[1] 3.50 0.62 2.45 4.45 1005 1.01
## alpha_var[2] 0.17 0.62 -0.86 1.12 996 1.01
## alpha_var[3] -2.85 0.84 -4.15 -1.48 1766 1.00
## alpha_var[4] 1.02 0.62 0.02 2.01 1005 1.01
## alpha_var[5] 1.20 0.61 0.21 2.19 981 1.01
## alpha_var[6] 1.10 0.61 0.15 2.10 962 1.01
## alpha_var[7] -0.68 0.64 -1.70 0.34 1062 1.01
## alpha_var[8] 0.46 0.61 -0.54 1.42 976 1.01
## alpha_var[9] 0.14 0.61 -0.87 1.10 981 1.01
## alpha_var[10] -3.24 0.85 -4.59 -1.90 1740 1.00
## alpha_var[11] 0.71 0.62 -0.26 1.72 984 1.01
## alpha_var[12] -0.09 0.62 -1.08 0.91 998 1.01
## alpha_var[13] 1.21 0.61 0.23 2.19 971 1.01
## alpha_var[14] -0.02 0.62 -1.00 0.98 993 1.01
## alpha_var[15] 1.75 0.61 0.76 2.71 966 1.01
## alpha_var[16] -2.54 0.80 -3.78 -1.24 1602 1.00
## alpha_var[17] -0.22 0.62 -1.26 0.73 1005 1.01
## alpha_var[18] -1.13 0.65 -2.18 -0.08 1120 1.01
## alpha_var[19] 0.14 0.62 -0.85 1.14 1003 1.01
## alpha_var[20] -1.90 0.67 -2.94 -0.79 1152 1.01
## alpha_var[21] -0.63 0.63 -1.62 0.38 1014 1.01
## alpha_var[22] -1.01 0.67 -2.07 0.08 1141 1.01
## alpha_var[23] 3.25 0.62 2.24 4.23 1001 1.01
## alpha_var[24] 0.00 0.62 -0.98 1.01 1006 1.01
## alpha_var[25] -4.72 1.20 -6.55 -2.84 2733 1.00
## alpha_var[26] -1.08 0.63 -2.08 -0.08 1009 1.01
## alpha_var[27] 0.41 0.62 -0.65 1.35 1005 1.01
## alpha_var[28] -0.02 0.62 -1.02 0.95 991 1.01
## alpha_var[29] -2.02 0.70 -3.12 -0.88 1290 1.00
## alpha_var[30] -1.54 0.66 -2.61 -0.49 1140 1.01
## alpha_var[31] 0.98 0.62 -0.03 1.96 997 1.01
## alpha_var[32] 0.03 0.62 -0.97 1.00 981 1.01
## alpha_var[33] -0.35 0.62 -1.35 0.66 1006 1.01
## alpha_var[34] 0.46 0.62 -0.55 1.42 991 1.01
## alpha_var[35] 1.52 0.61 0.54 2.49 969 1.01
## alpha_var[36] -1.22 0.68 -2.30 -0.13 1185 1.01
## alpha_var[37] -0.27 0.62 -1.30 0.69 997 1.01
## alpha_var[38] -1.23 0.65 -2.34 -0.25 1109 1.01
## alpha_var[39] -2.16 0.68 -3.23 -1.05 1207 1.01
## alpha_var[40] 1.32 0.61 0.35 2.31 975 1.01
## alpha_var[41] -1.82 0.66 -2.90 -0.77 1108 1.01
## alpha_var[42] 0.84 0.61 -0.15 1.82 978 1.01
## alpha_var[43] -0.22 0.63 -1.23 0.78 1015 1.01
## alpha_var[44] 0.97 0.62 0.00 1.96 990 1.01
## alpha_var[45] 1.44 0.62 0.46 2.43 982 1.01compare(m.cane.alpha,m.var)## WAIC pWAIC dWAIC weight SE dSE
## m.var 15758.6 45.5 0.0 1 170.53 NA
## m.cane.alpha 19489.1 1.0 3730.5 0 173.25 115.49
m.var.block <- map2stan(
alist(
r ~ dbinom(n,p),
logit(p) <- alpha + alpha_var[var_id] + alpha_block[block_id],
alpha <- dnorm(-1.55,1),
alpha_var[var_id] ~ dnorm(0,5),
alpha_block[block_id] ~ dnorm(0,sigma_block),
sigma_block ~ dcauchy(0,1)
),
data=cane,
chains = 4,
cores = 2,
iter=10000,
warmup=4000
)## Warning in FUN(X[[i]], ...): data with name var is not numeric and not used## Warning in FUN(X[[i]], ...): data with name block is not numeric and not
## used## Warning in FUN(X[[i]], ...): data with name var is not numeric and not used## Warning in FUN(X[[i]], ...): data with name block is not numeric and not
## used##
## SAMPLING FOR MODEL 'r ~ dbinom(n, p)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 4e-06 seconds (Warm-up)
## 0.000176 seconds (Sampling)
## 0.00018 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 2400 / 24000 ]
[ 4800 / 24000 ]
[ 7200 / 24000 ]
[ 9600 / 24000 ]
[ 12000 / 24000 ]
[ 14400 / 24000 ]
[ 16800 / 24000 ]
[ 19200 / 24000 ]
[ 21600 / 24000 ]
[ 24000 / 24000 ]## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.save.image("tomato1212.Rdata")